home *** CD-ROM | disk | FTP | other *** search
/ Workbench Add-On / Workbench Add-On - Volume 1.iso / Gfx / Edit / TSMorph / src / jpeg_ls / jrevdct.c < prev    next >
C/C++ Source or Header  |  1994-10-30  |  14KB  |  370 lines
  1. /*
  2.  * jrevdct.c
  3.  *
  4.  * Copyright (C) 1991, 1992, Thomas G. Lane.
  5.  * This file is part of the Independent JPEG Group's software.
  6.  * For conditions of distribution and use, see the accompanying README file.
  7.  *
  8.  * This file contains the basic inverse-DCT transformation subroutine.
  9.  *
  10.  * This implementation is based on an algorithm described in
  11.  *   C. Loeffler, A. Ligtenberg and G. Moschytz, "Practical Fast 1-D DCT
  12.  *   Algorithms with 11 Multiplications", Proc. Int'l. Conf. on Acoustics,
  13.  *   Speech, and Signal Processing 1989 (ICASSP '89), pp. 988-991.
  14.  * The primary algorithm described there uses 11 multiplies and 29 adds.
  15.  * We use their alternate method with 12 multiplies and 32 adds.
  16.  * The advantage of this method is that no data path contains more than one
  17.  * multiplication; this allows a very simple and accurate implementation in
  18.  * scaled fixed-point arithmetic, with a minimal number of shifts.
  19.  */
  20.  
  21. #include "jinclude.h"
  22.  
  23. /*
  24.  * This routine is specialized to the case DCTSIZE = 8.
  25.  */
  26.  
  27. #if DCTSIZE != 8
  28.   Sorry, this code only copes with 8x8 DCTs. /* deliberate syntax err */
  29. #endif
  30.  
  31.  
  32. /*
  33.  * A 2-D IDCT can be done by 1-D IDCT on each row followed by 1-D IDCT
  34.  * on each column.  Direct algorithms are also available, but they are
  35.  * much more complex and seem not to be any faster when reduced to code.
  36.  *
  37.  * The poop on this scaling stuff is as follows:
  38.  *
  39.  * Each 1-D IDCT step produces outputs which are a factor of sqrt(N)
  40.  * larger than the true IDCT outputs.  The final outputs are therefore
  41.  * a factor of N larger than desired; since N=8 this can be cured by
  42.  * a simple right shift at the end of the algorithm.  The advantage of
  43.  * this arrangement is that we save two multiplications per 1-D IDCT,
  44.  * because the y0 and y4 inputs need not be divided by sqrt(N).
  45.  *
  46.  * We have to do addition and subtraction of the integer inputs, which
  47.  * is no problem, and multiplication by fractional constants, which is
  48.  * a problem to do in integer arithmetic.  We multiply all the constants
  49.  * by CONST_SCALE and convert them to integer constants (thus retaining
  50.  * CONST_BITS bits of precision in the constants).  After doing a
  51.  * multiplication we have to divide the product by CONST_SCALE, with proper
  52.  * rounding, to produce the correct output.  This division can be done
  53.  * cheaply as a right shift of CONST_BITS bits.  We postpone shifting
  54.  * as long as possible so that partial sums can be added together with
  55.  * full fractional precision.
  56.  *
  57.  * The outputs of the first pass are scaled up by PASS1_BITS bits so that
  58.  * they are represented to better-than-integral precision.  These outputs
  59.  * require BITS_IN_JSAMPLE + PASS1_BITS + 3 bits; this fits in a 16-bit word
  60.  * with the recommended scaling.  (To scale up 12-bit sample data further, an
  61.  * intermediate INT32 array would be needed.)
  62.  *
  63.  * To avoid overflow of the 32-bit intermediate results in pass 2, we must
  64.  * have BITS_IN_JSAMPLE + CONST_BITS + PASS1_BITS <= 26.  Error analysis
  65.  * shows that the values given below are the most effective.
  66.  */
  67.  
  68. #ifdef EIGHT_BIT_SAMPLES
  69. #define CONST_BITS  13
  70. #define PASS1_BITS  2
  71. #else
  72. #define CONST_BITS  13
  73. #define PASS1_BITS  1        /* lose a little precision to avoid overflow */
  74. #endif
  75.  
  76. #define ONE    ((INT32) 1)
  77.  
  78. #define CONST_SCALE (ONE << CONST_BITS)
  79.  
  80. /* Convert a positive real constant to an integer scaled by CONST_SCALE. */
  81.  
  82. #define FIX(x)    ((INT32) ((x) * CONST_SCALE + 0.5))
  83.  
  84. /* Some C compilers fail to reduce "FIX(constant)" at compile time, thus
  85.  * causing a lot of useless floating-point operations at run time.
  86.  * To get around this we use the following pre-calculated constants.
  87.  * If you change CONST_BITS you may want to add appropriate values.
  88.  * (With a reasonable C compiler, you can just rely on the FIX() macro...)
  89.  */
  90.  
  91. #if CONST_BITS == 13
  92. #define FIX_0_298631336  ((INT32)  2446)    /* FIX(0.298631336) */
  93. #define FIX_0_390180644  ((INT32)  3196)    /* FIX(0.390180644) */
  94. #define FIX_0_541196100  ((INT32)  4433)    /* FIX(0.541196100) */
  95. #define FIX_0_765366865  ((INT32)  6270)    /* FIX(0.765366865) */
  96. #define FIX_0_899976223  ((INT32)  7373)    /* FIX(0.899976223) */
  97. #define FIX_1_175875602  ((INT32)  9633)    /* FIX(1.175875602) */
  98. #define FIX_1_501321110  ((INT32)  12299)    /* FIX(1.501321110) */
  99. #define FIX_1_847759065  ((INT32)  15137)    /* FIX(1.847759065) */
  100. #define FIX_1_961570560  ((INT32)  16069)    /* FIX(1.961570560) */
  101. #define FIX_2_053119869  ((INT32)  16819)    /* FIX(2.053119869) */
  102. #define FIX_2_562915447  ((INT32)  20995)    /* FIX(2.562915447) */
  103. #define FIX_3_072711026  ((INT32)  25172)    /* FIX(3.072711026) */
  104. #else
  105. #define FIX_0_298631336  FIX(0.298631336)
  106. #define FIX_0_390180644  FIX(0.390180644)
  107. #define FIX_0_541196100  FIX(0.541196100)
  108. #define FIX_0_765366865  FIX(0.765366865)
  109. #define FIX_0_899976223  FIX(0.899976223)
  110. #define FIX_1_175875602  FIX(1.175875602)
  111. #define FIX_1_501321110  FIX(1.501321110)
  112. #define FIX_1_847759065  FIX(1.847759065)
  113. #define FIX_1_961570560  FIX(1.961570560)
  114. #define FIX_2_053119869  FIX(2.053119869)
  115. #define FIX_2_562915447  FIX(2.562915447)
  116. #define FIX_3_072711026  FIX(3.072711026)
  117. #endif
  118.  
  119.  
  120. /* Descale and correctly round an INT32 value that's scaled by N bits.
  121.  * We assume RIGHT_SHIFT rounds towards minus infinity, so adding
  122.  * the fudge factor is correct for either sign of X.
  123.  */
  124.  
  125. #define DESCALE(x,n)  RIGHT_SHIFT((x) + (ONE << ((n)-1)), n)
  126.  
  127. /* Multiply an INT32 variable by an INT32 constant to yield an INT32 result.
  128.  * For 8-bit samples with the recommended scaling, all the variable
  129.  * and constant values involved are no more than 16 bits wide, so a
  130.  * 16x16->32 bit multiply can be used instead of a full 32x32 multiply;
  131.  * this provides a useful speedup on many machines.
  132.  * There is no way to specify a 16x16->32 multiply in portable C, but
  133.  * some C compilers will do the right thing if you provide the correct
  134.  * combination of casts.
  135.  * NB: for 12-bit samples, a full 32-bit multiplication will be needed.
  136.  */
  137.  
  138. #ifdef EIGHT_BIT_SAMPLES
  139. #ifdef SHORTxSHORT_32        /* may work if 'int' is 32 bits */
  140. #define MULTIPLY(var,const)  (((INT16) (var)) * ((INT16) (const)))
  141. #endif
  142. #ifdef SHORTxLCONST_32        /* known to work with Microsoft C 6.0 */
  143. #define MULTIPLY(var,const)  (((INT16) (var)) * ((INT32) (const)))
  144. #endif
  145. #endif
  146.  
  147. #ifndef MULTIPLY        /* default definition */
  148. #define MULTIPLY(var,const)  ((var) * (const))
  149. #endif
  150.  
  151.  
  152. /*
  153.  * Perform the inverse DCT on one block of coefficients.
  154.  */
  155.  
  156. GLOBAL void
  157. j_rev_dct (DCTBLOCK data)
  158. {
  159.   INT32 tmp0, tmp1, tmp2, tmp3;
  160.   INT32 tmp10, tmp11, tmp12, tmp13;
  161.   INT32 z1, z2, z3, z4, z5;
  162.   register DCTELEM *dataptr;
  163.   int rowctr;
  164.   SHIFT_TEMPS
  165.  
  166.   /* Pass 1: process rows. */
  167.   /* Note results are scaled up by sqrt(8) compared to a true IDCT; */
  168.   /* furthermore, we scale the results by 2**PASS1_BITS. */
  169.  
  170.   dataptr = data;
  171.   for (rowctr = DCTSIZE-1; rowctr >= 0; rowctr--) {
  172.     /* Due to quantization, we will usually find that many of the input
  173.      * coefficients are zero, especially the AC terms.  We can exploit this
  174.      * by short-circuiting the IDCT calculation for any row in which all
  175.      * the AC terms are zero.  In that case each output is equal to the
  176.      * DC coefficient (with scale factor as needed).
  177.      * With typical images and quantization tables, half or more of the
  178.      * row DCT calculations can be simplified this way.
  179.      */
  180.  
  181.     if ((dataptr[1] | dataptr[2] | dataptr[3] | dataptr[4] |
  182.      dataptr[5] | dataptr[6] | dataptr[7]) == 0) {
  183.       /* AC terms all zero */
  184.       DCTELEM dcval = (DCTELEM) (dataptr[0] << PASS1_BITS);
  185.       
  186.       dataptr[0] = dcval;
  187.       dataptr[1] = dcval;
  188.       dataptr[2] = dcval;
  189.       dataptr[3] = dcval;
  190.       dataptr[4] = dcval;
  191.       dataptr[5] = dcval;
  192.       dataptr[6] = dcval;
  193.       dataptr[7] = dcval;
  194.       
  195.       dataptr += DCTSIZE;    /* advance pointer to next row */
  196.       continue;
  197.     }
  198.  
  199.     /* Even part: reverse the even part of the forward DCT. */
  200.     /* The rotator is sqrt(2)*c(-6). */
  201.  
  202.     z2 = (INT32) dataptr[2];
  203.     z3 = (INT32) dataptr[6];
  204.  
  205.     z1 = MULTIPLY(z2 + z3, FIX_0_541196100);
  206.     tmp2 = z1 + MULTIPLY(z3, - FIX_1_847759065);
  207.     tmp3 = z1 + MULTIPLY(z2, FIX_0_765366865);
  208.  
  209.     tmp0 = ((INT32) dataptr[0] + (INT32) dataptr[4]) << CONST_BITS;
  210.     tmp1 = ((INT32) dataptr[0] - (INT32) dataptr[4]) << CONST_BITS;
  211.  
  212.     tmp10 = tmp0 + tmp3;
  213.     tmp13 = tmp0 - tmp3;
  214.     tmp11 = tmp1 + tmp2;
  215.     tmp12 = tmp1 - tmp2;
  216.     
  217.     /* Odd part per figure 8; the matrix is unitary and hence its
  218.      * transpose is its inverse.  i0..i3 are y7,y5,y3,y1 respectively.
  219.      */
  220.  
  221.     tmp0 = (INT32) dataptr[7];
  222.     tmp1 = (INT32) dataptr[5];
  223.     tmp2 = (INT32) dataptr[3];
  224.     tmp3 = (INT32) dataptr[1];
  225.  
  226.     z1 = tmp0 + tmp3;
  227.     z2 = tmp1 + tmp2;
  228.     z3 = tmp0 + tmp2;
  229.     z4 = tmp1 + tmp3;
  230.     z5 = MULTIPLY(z3 + z4, FIX_1_175875602); /* sqrt(2) * c3 */
  231.     
  232.     tmp0 = MULTIPLY(tmp0, FIX_0_298631336); /* sqrt(2) * (-c1+c3+c5-c7) */
  233.     tmp1 = MULTIPLY(tmp1, FIX_2_053119869); /* sqrt(2) * ( c1+c3-c5+c7) */
  234.     tmp2 = MULTIPLY(tmp2, FIX_3_072711026); /* sqrt(2) * ( c1+c3+c5-c7) */
  235.     tmp3 = MULTIPLY(tmp3, FIX_1_501321110); /* sqrt(2) * ( c1+c3-c5-c7) */
  236.     z1 = MULTIPLY(z1, - FIX_0_899976223); /* sqrt(2) * (c7-c3) */
  237.     z2 = MULTIPLY(z2, - FIX_2_562915447); /* sqrt(2) * (-c1-c3) */
  238.     z3 = MULTIPLY(z3, - FIX_1_961570560); /* sqrt(2) * (-c3-c5) */
  239.     z4 = MULTIPLY(z4, - FIX_0_390180644); /* sqrt(2) * (c5-c3) */
  240.     
  241.     z3 += z5;
  242.     z4 += z5;
  243.     
  244.     tmp0 += z1 + z3;
  245.     tmp1 += z2 + z4;
  246.     tmp2 += z2 + z3;
  247.     tmp3 += z1 + z4;
  248.  
  249.     /* Final output stage: inputs are tmp10..tmp13, tmp0..tmp3 */
  250.  
  251.     dataptr[0] = (DCTELEM) DESCALE(tmp10 + tmp3, CONST_BITS-PASS1_BITS);
  252.     dataptr[7] = (DCTELEM) DESCALE(tmp10 - tmp3, CONST_BITS-PASS1_BITS);
  253.     dataptr[1] = (DCTELEM) DESCALE(tmp11 + tmp2, CONST_BITS-PASS1_BITS);
  254.     dataptr[6] = (DCTELEM) DESCALE(tmp11 - tmp2, CONST_BITS-PASS1_BITS);
  255.     dataptr[2] = (DCTELEM) DESCALE(tmp12 + tmp1, CONST_BITS-PASS1_BITS);
  256.     dataptr[5] = (DCTELEM) DESCALE(tmp12 - tmp1, CONST_BITS-PASS1_BITS);
  257.     dataptr[3] = (DCTELEM) DESCALE(tmp13 + tmp0, CONST_BITS-PASS1_BITS);
  258.     dataptr[4] = (DCTELEM) DESCALE(tmp13 - tmp0, CONST_BITS-PASS1_BITS);
  259.  
  260.     dataptr += DCTSIZE;        /* advance pointer to next row */
  261.   }
  262.  
  263.   /* Pass 2: process columns. */
  264.   /* Note that we must descale the results by a factor of 8 == 2**3, */
  265.   /* and also undo the PASS1_BITS scaling. */
  266.  
  267.   dataptr = data;
  268.   for (rowctr = DCTSIZE-1; rowctr >= 0; rowctr--) {
  269.     /* Columns of zeroes can be exploited in the same way as we did with rows.
  270.      * However, the row calculation has created many nonzero AC terms, so the
  271.      * simplification applies less often (typically 5% to 10% of the time).
  272.      * On machines with very fast multiplication, it's possible that the
  273.      * test takes more time than it's worth.  In that case this section
  274.      * may be commented out.
  275.      */
  276.  
  277. #ifndef NO_ZERO_COLUMN_TEST
  278.     if ((dataptr[DCTSIZE*1] | dataptr[DCTSIZE*2] | dataptr[DCTSIZE*3] |
  279.      dataptr[DCTSIZE*4] | dataptr[DCTSIZE*5] | dataptr[DCTSIZE*6] |
  280.      dataptr[DCTSIZE*7]) == 0) {
  281.       /* AC terms all zero */
  282.       DCTELEM dcval = (DCTELEM) DESCALE((INT32) dataptr[0], PASS1_BITS+3);
  283.       
  284.       dataptr[DCTSIZE*0] = dcval;
  285.       dataptr[DCTSIZE*1] = dcval;
  286.       dataptr[DCTSIZE*2] = dcval;
  287.       dataptr[DCTSIZE*3] = dcval;
  288.       dataptr[DCTSIZE*4] = dcval;
  289.       dataptr[DCTSIZE*5] = dcval;
  290.       dataptr[DCTSIZE*6] = dcval;
  291.       dataptr[DCTSIZE*7] = dcval;
  292.       
  293.       dataptr++;        /* advance pointer to next column */
  294.       continue;
  295.     }
  296. #endif
  297.  
  298.     /* Even part: reverse the even part of the forward DCT. */
  299.     /* The rotator is sqrt(2)*c(-6). */
  300.  
  301.     z2 = (INT32) dataptr[DCTSIZE*2];
  302.     z3 = (INT32) dataptr[DCTSIZE*6];
  303.  
  304.     z1 = MULTIPLY(z2 + z3, FIX_0_541196100);
  305.     tmp2 = z1 + MULTIPLY(z3, - FIX_1_847759065);
  306.     tmp3 = z1 + MULTIPLY(z2, FIX_0_765366865);
  307.  
  308.     tmp0 = ((INT32) dataptr[DCTSIZE*0] + (INT32) dataptr[DCTSIZE*4]) << CONST_BITS;
  309.     tmp1 = ((INT32) dataptr[DCTSIZE*0] - (INT32) dataptr[DCTSIZE*4]) << CONST_BITS;
  310.  
  311.     tmp10 = tmp0 + tmp3;
  312.     tmp13 = tmp0 - tmp3;
  313.     tmp11 = tmp1 + tmp2;
  314.     tmp12 = tmp1 - tmp2;
  315.     
  316.     /* Odd part per figure 8; the matrix is unitary and hence its
  317.      * transpose is its inverse.  i0..i3 are y7,y5,y3,y1 respectively.
  318.      */
  319.  
  320.     tmp0 = (INT32) dataptr[DCTSIZE*7];
  321.     tmp1 = (INT32) dataptr[DCTSIZE*5];
  322.     tmp2 = (INT32) dataptr[DCTSIZE*3];
  323.     tmp3 = (INT32) dataptr[DCTSIZE*1];
  324.  
  325.     z1 = tmp0 + tmp3;
  326.     z2 = tmp1 + tmp2;
  327.     z3 = tmp0 + tmp2;
  328.     z4 = tmp1 + tmp3;
  329.     z5 = MULTIPLY(z3 + z4, FIX_1_175875602); /* sqrt(2) * c3 */
  330.     
  331.     tmp0 = MULTIPLY(tmp0, FIX_0_298631336); /* sqrt(2) * (-c1+c3+c5-c7) */
  332.     tmp1 = MULTIPLY(tmp1, FIX_2_053119869); /* sqrt(2) * ( c1+c3-c5+c7) */
  333.     tmp2 = MULTIPLY(tmp2, FIX_3_072711026); /* sqrt(2) * ( c1+c3+c5-c7) */
  334.     tmp3 = MULTIPLY(tmp3, FIX_1_501321110); /* sqrt(2) * ( c1+c3-c5-c7) */
  335.     z1 = MULTIPLY(z1, - FIX_0_899976223); /* sqrt(2) * (c7-c3) */
  336.     z2 = MULTIPLY(z2, - FIX_2_562915447); /* sqrt(2) * (-c1-c3) */
  337.     z3 = MULTIPLY(z3, - FIX_1_961570560); /* sqrt(2) * (-c3-c5) */
  338.     z4 = MULTIPLY(z4, - FIX_0_390180644); /* sqrt(2) * (c5-c3) */
  339.     
  340.     z3 += z5;
  341.     z4 += z5;
  342.     
  343.     tmp0 += z1 + z3;
  344.     tmp1 += z2 + z4;
  345.     tmp2 += z2 + z3;
  346.     tmp3 += z1 + z4;
  347.  
  348.     /* Final output stage: inputs are tmp10..tmp13, tmp0..tmp3 */
  349.  
  350.     dataptr[DCTSIZE*0] = (DCTELEM) DESCALE(tmp10 + tmp3,
  351.                        CONST_BITS+PASS1_BITS+3);
  352.     dataptr[DCTSIZE*7] = (DCTELEM) DESCALE(tmp10 - tmp3,
  353.                        CONST_BITS+PASS1_BITS+3);
  354.     dataptr[DCTSIZE*1] = (DCTELEM) DESCALE(tmp11 + tmp2,
  355.                        CONST_BITS+PASS1_BITS+3);
  356.     dataptr[DCTSIZE*6] = (DCTELEM) DESCALE(tmp11 - tmp2,
  357.                        CONST_BITS+PASS1_BITS+3);
  358.     dataptr[DCTSIZE*2] = (DCTELEM) DESCALE(tmp12 + tmp1,
  359.                        CONST_BITS+PASS1_BITS+3);
  360.     dataptr[DCTSIZE*5] = (DCTELEM) DESCALE(tmp12 - tmp1,
  361.                        CONST_BITS+PASS1_BITS+3);
  362.     dataptr[DCTSIZE*3] = (DCTELEM) DESCALE(tmp13 + tmp0,
  363.                        CONST_BITS+PASS1_BITS+3);
  364.     dataptr[DCTSIZE*4] = (DCTELEM) DESCALE(tmp13 - tmp0,
  365.                        CONST_BITS+PASS1_BITS+3);
  366.     
  367.     dataptr++;            /* advance pointer to next column */
  368.   }
  369. }
  370.